\documentclass{article}

\usepackage{graphicx} % Required for inserting images
\usepackage{amssymb}
\usepackage{mathrsfs}
\usepackage{amsmath}
\usepackage{mathrsfs}
\usepackage{xltxtra}
\usepackage{listings}
\usepackage{graphicx} % Required for inserting images
\usepackage{amsmath}
\usepackage{xeCJK}
\setCJKmainfont{Noto Serif CJK TC}

\title{Proj}
\author{王薪雅 3220300210}
\date{\today}

\begin{document}

\maketitle

\section{介绍}
\textbf{有限元方法的基本设置}\\
这是我们实际使用有限元来计算某些内容的第一个示例。 我们将求解边界值为零但右侧非零的泊松方程的简单版本：\\
$
\begin{array}{ccc}
    -\Delta u & = f & \text{in } \Omega, \\
    u & = 0 &   \text{on } \partial \Omega.
\end{array}
$
\\

我们将在正方形 $\Omega= [−1,1]^2$上求解这个方程，您已经了解了如何在步骤1和步骤2中生成网格。在此程序中，我们也将仅考虑特定情况 $f(x)=1$。

如果您已经了解了有限元方法的基础知识，您将记住我们通过有限维近似来近似解 $u$ 所需的步骤。 具体来说，我们首先需要推导出上面方程的弱形式，我们通过将方程乘以左侧的测试函数 $\varphi$ 来获得（下面我们会回到从左侧相乘而不是从右侧相乘的原因） 并在域 $\Omega$ 上积分：\\
$-\int_{\Omega} \varphi \Delta u=\int_{\Omega} \varphi f$ \\
这可以通过部分积分：\\
$\int_{\Omega} \nabla \varphi \cdot \nabla u-\int{\partial \Omega}n \cdot \nabla u = \int_{\Omega}\varphi f$.\\
测试函数$\varphi$必须满足相同类型的边界条件（用数学术语来说：它需要来自我们寻求解的集合的切空间），因此在边界 $\varphi=0$ 上，因此弱形式我们 正在寻找读物$(\nabla \varphi,\nabla u)=(\varphi,f)$,
其中我们使用了通用符号 $(a,b)=\int_{\Omega}a b$。 然后，问题要求一个函数 $u$，对于来自适当空间（这里是空间 $H^1$）的所有测试函数$\varphi$，该陈述都成立。

当然，一般情况下我们在计算机上找不到这样的函数，而是寻求一个近似值 $u_h(x)=\Sigma_j U_j\varphi_j(x)$，其中 $U_j$ 是我们需要确定的未知展开系数（“度”） 这个问题的自由度”），而 $\varphi_i(x)$ 是我们将使用的有限元形函数。 为了定义这些形状函数，我们需要以下内容：
\begin{itemize}
    \item 用于定义形状函数的网格。 您已经在步骤 1 和步骤 2 中了解了如何生成和操作描述网格的对象。
    \item 描述我们想要在参考单元上使用的形状函数的有限元（在 deal.II 中始终是单位间隔 [0,1]、单位平方 $[0,1]^2$ 或单位立方体$[0,1]^3$，取决于您工作的空间维度）。 在步骤 2 中，我们已经使用了 FE\_Q<2>类型的对象，它表示通过在支撑点上插值来定义形状函数的常见拉格朗日元素。 最简单的是FE\_Q<2>(1)，它使用多项式次数1。在2d中，这些通常被称为双线性，因为它们在参考单元的两个坐标中的每一个中都是线性的。 （在 1d 中，它们是线性的，在 3d 中是三线性的；但是，在 deal.II 文档中，我们通常不会进行这种区分，而只是始终将这些函数称为“线性”。）
    \item 一个 DoFHandler 对象，以有限元对象提供的参考单元描述为基础，枚举网格上的所有自由度。 您也已经在步骤 2 中了解了如何执行此操作。
    \item 一种映射，说明如何从参考单元上的有限元类定义的形状函数获得实际单元上的形状函数。 默认情况下，除非您明确另有说明，否则 deal.II 将为此使用（双、三）线性映射，因此在大多数情况下您不必担心此步骤。
\end{itemize}
通过这些步骤，我们现在有了一组函数 $\varphi_i$，我们可以定义离散问题的弱形式： 找到一个函数 $u_h$，即找到上面提到的展开系数 $U_j$,使得\\
$(\nabla \varphi_i,\nabla u_h)=(\varphi_i,f)$  ,$i=0\dots N-1$.\\
请注意，我们这里遵循的约定是所有内容都从零开始计数，这在 C 和 C++ 中很常见。 如果插入表示 $u_h(x)=\Sigma_j U_j \varphi_j(x)$ 然后观察到，该方程可以重写为线性系统\\
$
\begin{aligned}
\left(\nabla \varphi_i, \nabla u_h\right) & =\left(\nabla \varphi_i, \nabla\left[\sum_j U_j \varphi_j\right]\right) \\
& =\sum_j\left(\nabla \varphi_i, \nabla\left[U_j \varphi_j\right]\right) \\
& =\sum_j\left(\nabla \varphi_i, \nabla \varphi_j\right) U_j .
\end{aligned}
$ \\
这样，问题就变成了：找到一个向量 $U$，使得$AU=F$,其中矩阵 $A$ 和右侧 $F$ 定义为\\
$
\begin{array}{cc}
    A_{ij} & =(\nabla \varphi_i,\nabla \varphi_j) \\
    F_i & =(\varphi_i,f).
\end{array}
$ \\
\textbf{我们应该从左边还是从右边乘以测试函数？}

在我们继续描述如何计算这些量之前，请注意，如果我们将原始方程从右侧而不是从左侧乘以测试函数，那么我们将获得以下形式的线性系统
$U^TA=F^T$
与行向量 $F^T$。 通过转置这个系统，这当然等价于求解$A^TU=F$。这里与上面相同，因为 $A=A^T$。 但总的来说并非如此，为了避免任何形式的混乱，经验表明，简单地养成从左而不是从右相乘方程的习惯（就像数学文献中经常做的那样）可以避免常见的情况。 错误类别，因为矩阵会自动正确，并且在比较理论和实现时不需要转置。 请参阅本教程中第一个示例的步骤9，其中我们有一个非对称双线性形式，无论从右侧还是从左侧相乘，都会有所不同。\\
\textbf{组装矩阵和右侧向量}

现在我们知道我们需要什么（即：保存矩阵和向量的对象，以及计算 $A_{ij},F_i$ 的方法），我们可以看看如何实现这一点：
\begin{itemize}
    \item $A$ 的对象属于 SparseMatrix 类型，而 $U$ 和 $F$ 的对象属于 Vector 类型。 我们将在下面的程序中看到哪些类用于求解线性系统。
    \item 我们需要一种形成积分的方法。在有限元方法中，最常见的是使用求积法来完成，即用每个单元上的一组求积点的加权和来代替积分。也就是说，我们首先将$\Omega$上的积分拆分为所有单元上的积分，\\
    $
    \begin{array}{lll}
         A_{ij} & = (\nabla \varphi_i,\nabla \varphi_j) & = \sum_{K \in \mathbb{T}} \int_K \nabla \varphi_i \cdot \nabla \varphi_j ,  \\
         F_i & = (\varphi_i,f) & = \sum_{K \in \mathbb{T}} \int_K \varphi_i f,
    \end{array}
    $

然后通过求积来近似每个单元格的贡献：\\
$
\begin{array}{ll}
     A_{ij}^K & = \int_K \nabla \varphi_i \cdot \nabla \varphi_j \approx \sum_q \nabla \varphi_i(x_q^K)\cdot \nabla \varphi_j(x_q^K)w_q^K,  \\
     F_i^K & = \int_K \varphi_i f\approx \sum_q \varphi_i(x_q^K)f(x_q^K)w_q^K, 
\end{array}
$

其中 $\mathbb{T} \approx \Omega$ 是近似域的三角剖分，$x^K_q$ 是单元 $K$ 上的第 $q$ 个正交点，$w^K_q$ 是第 $q$ 个正交权重。 执行此操作需要不同的部分，我们将在接下来依次讨论它们。
\item 首先，我们需要一种方法来描述正交点的位置$x^K_q$及其权重$w^K_q$。 它们通常以与形状函数相同的方式从参考单元映射，即隐式使用 MappingQ1 类，或者，如果您明确这么说，则通过从 Mapping 派生的其他类之一。 参考单元上的位置和权重由派生自 Quadrature 基类的对象描述。 通常，选择求积公式（即一组点和权重），使求积完全等于矩阵中的积分； 这是可以实现的，因为积分中的所有因子都是多项式，并且是通过高斯求积公式完成的，在 QGauss 类中实现。
\item 然后，我们需要一些可以帮助我们计算单元 $K$ 上的 $\varphi_i(x_q^K)$ 的东西。这就是 FEValues 类的作用：它需要一个有限元对象来描述参考单元上的$\varphi$，一个求积对象来描述正交点和权重， 和一个映射对象（或隐式采用 MappingQ1 类），并提供真实单元 $K$ 上形状函数的值和导数，以及位于 $K$ 上的正交点处积分所需的各种其他信息。
\end{itemize}

将矩阵和右侧计算为所有单元格的总和（然后计算正交点的总和）的过程通常称为组装线性系统，或简称组装，使用与装配线相关的单词的含义，意思是 “将一组碎片、片段或元素组合在一起的行为”。

FEValues 确实是组装过程中的中心类。 一种查看它的方法如下：FiniteElement 和派生类描述形状函数，即无限维对象：函数在每个点都有值。 出于理论原因，我们需要这个，因为我们想要使用函数积分来执行分析。 然而，对于计算机来说，这是一个非常困难的概念，因为它们通常只能处理有限量的信息，因此我们用通过使用定义的点进行映射（Mapping 对象）获得的积分点的和来替换积分 将参考单元（求积对象）上的点投影到真实单元上的点上。 本质上，我们将问题简化为只需要有限数量的信息，即形状函数值和导数、正交权重、法向量等，仅在有限的点集上。 FEValues 类将三个组件组合在一起，并提供有关特定单元 $K$ 的有限信息集。当我们组装下面的线性系统时，您将看到它的实际作用。

值得注意的是，如果您只是在应用程序中自己创建这三个对象，并自己处理信息，那么所有这一切也可以实现。 然而，这既不会更简单（FEValues 类提供了您实际需要的信息）也不会更快：FEValues 类经过高度优化，仅在每个单元上计算您需要的特定信息； 如果可以从前一个单元中重用任何内容，那么它将这样做，并且该类中有很多代码来确保在有利的地方缓存内容。

本介绍的最后一部分是要提到，在获得线性系统后，使用迭代求解器对其进行求解，然后进行后处理：我们使用 DataOut 类创建一个输出文件，然后可以使用一种常见的可视化程序对其进行可视化。

\textbf{求解线性系统}

对于有限元程序，我们最终得到的线性系统相对较小：矩阵的大小为 $1089\times1089$，因为我们使用的网格是 $32\times32$，因此网格中有 $332=1089$ 个顶点 。 在后来的许多教程程序中，矩阵大小在数万到数十万的范围内并不少见，并且使用诸如ASPECT之类基于deal.II构建的代码，我们经常解决超过一亿个方程的问题 （尽管使用并行计算机）。 无论如何，即使对于这里的小系统，矩阵也比本科生或大多数研究生课程中通常遇到的矩阵大得多，因此出现了我们如何解决此类线性系统的问题。

人们通常学习的用于求解线性系统的第一种方法是高斯消去法。 这种方法的问题在于，它需要与 $N^3$ 成正比的运算次数，其中 $N$ 是线性系统中方程或未知数的数量 - 更具体地说，运算次数为 $\frac{2}{3}N^3$，或多或少。 当 $N=1089$ 时，这意味着我们必须执行大约 8.61 亿次操作。 这是一个相当可行的数字，现代处理器只需不到 0.1 秒即可完成此操作。 但很明显，这不会扩展：如果线性系统中的方程数量是原来的 20 倍（即未知数的 20 倍），那么就已经需要 1000-10,000 秒或大约为 一小时。 将线性系统再放大十倍，显然我们无法再在一台计算机上解决它。

通过认识到矩阵中只有相对少量的条目是非零的，即矩阵是稀疏的，我们可以在某种程度上挽救这种情况。 高斯消去法的变体可以利用这一点，使过程大大加快； 我们将在步骤 29 中首次使用一个这样的方法（在 SparseDirectUMFPACK 类中实现），以及此后的其他方法。 高斯消元法的这些变化可能会让我们的问题规模达到 100,000 或 200,000 的量级，但也不会超出这个范围。

相反，我们在这里要做的是采用 1952 年的一个想法：共轭梯度法，简称“CG”。 CG 是一个“迭代”求解器，因为它形成了收敛到精确解的向量序列； 事实上，在没有舍入误差的情况下，经过 $N$ 次这样的迭代后，如果矩阵是对称且正定的，它会找到精确解。 该方法最初是作为另一种精确求解线性系统的方法而开发的，就像高斯消元法一样，但因此它几乎没有什么优点，并且在几十年内基本上被遗忘了。 但是，当计算机变得强大到足以解决高斯消去法不再适用的问题时（在 20 世纪 80 年代的某个时候），CG 被重新发现，因为人们意识到它非常适合像我们这样的大型稀疏系统。 由有限元法得到。 这是因为 (i) 它计算的向量收敛于精确解，因此只要我们对相当好的近似值感到满意，我们实际上不必进行所有 $N$ 次迭代来找到精确解； (ii) 它只需要矩阵向量乘积，这对于稀疏矩阵非常有用，因为根据定义，稀疏矩阵只有 $\mathscr{O}(N)$ 个条目，因此可以使用 $\mathscr{O}(N)$ 完成矩阵向量乘积 需要花费 $N^2$ 运算才能对稠密矩阵执行相同操作。 因此，我们希望用最多 $\mathscr{O}(N^2)$ 次运算来求解线性系统，并且在许多情况下要少得多。

因此，有限元代码几乎总是使用迭代求解器（例如 CG）来求解线性系统，我们也将在此代码中这样做。 （我们注意到，CG 方法仅适用于对称和正定矩阵；对于其他方程，矩阵可能不具有这些属性，我们将不得不使用迭代求解器的其他变体，例如适用于 更一般的矩阵。）

这些迭代求解器的一个重要组成部分是我们指定我们想要求解线性系统的容差——本质上，这是关于我们愿意在近似解中接受的误差的声明。 线性系统 $Ax=b$ 的精确解 $x$ 得到的近似解 $\Tilde{x}$ 的误差被定义为 $\Arrowvert x-\Tilde{x} \Arrowvert$，但这是一个我们无法计算的量，因为我们不知道精确解 $x$。 相反，我们通常将残差（定义为 $\Arrowvert b-A\Tilde{x} \Arrowvert = \Arrowvert A(x-\Tilde{x})\Arrowvert$视为可计算的度量。 然后我们让迭代求解器计算出越来越精确的解 $\Tilde{x}$，直到$\Arrowvert b-A\Tilde{x} \Arrowvert \leqslant \tau$ 。 一个实际问题是 τ 应该取什么值。 在大多数应用中，设置$\tau = 10^{-6}\Arrowvert b \Arrowvert$
是一个合理的选择。 事实上，我们使$\tau$与 $b$ 的大小（范数）成正比，这确保了我们对解的精度的预期与解的大小相关。 这是有道理的：如果我们将右侧 $b$ 放大十倍，那么 $Ax=b$ 的解 $x$ 也会变大十倍，$\Tilde{x}$ 也会变大十倍； 我们希望 $\Tilde{x}$ 中的精确数字数量与之前相同，这意味着当残差 $\Arrowvert b-A\Tilde{x}\Arrowvert$ 是原始大小的十倍时，我们也应该终止 - 这正是我们使 $\tau$ 与 $\Arrowvert b \Arrowvert$ 成比例时得到的结果。

所有这些都将在该程序的 Step3::solve() 函数中实现。 正如您将看到的，使用 deal.II 设置线性求解器非常简单：整个函数只有三行。

\textbf{关于实施}

尽管这是可以使用有限元方法求解的最简单的方程，但该程序显示了大多数有限元程序的基本结构，并且也充当了几乎所有以下程序都将遵循的模板。 具体来说，该程序的主类如下所示：
\begin{lstlisting}
class Step3
{
  public:
    Step3 ();
    void run ();
 
  private:
    void make_grid ();
    void setup_system ();
    void assemble_system ();
    void solve ();
    void output_results () const;
 
    Triangulation<2>     triangulation;
    FE_Q<2>              fe;
    DoFHandler<2>        dof_handler;
 
    SparsityPattern      sparsity_pattern;
    SparseMatrix<double> system_matrix;
    Vector<double>       solution;
    Vector<double>       system_rhs;
};
\end{lstlisting}

这遵循面向对象编程的数据封装原则，即我们尽力将此类的几乎所有内部细节隐藏在外部无法访问的私有成员中。

让我们从成员变量开始：这些变量遵循我们上面在要点中概述的构建块，即我们需要一个 Triangulation 和一个 DoFHandler 对象，以及一个描述我们想要使用的形状函数类型的有限元对象。 第二组对象与线性代数相关：系统矩阵和右侧以及解向量，以及描述矩阵稀疏模式的对象。 这就是此类的全部需求（以及任何平稳偏微分方程求解器所需的要素），并且需要在整个程序中生存。 与此相反，我们组装所需的 FEValues 对象仅在整个组装过程中需要，因此我们在执行此操作的函数中将其创建为本地对象，并在结束时再次销毁它。

其次，我们看一下成员函数。 这些也已经形成了几乎所有以下教程程序都会使用的通用结构：
\begin{itemize}
    \item make\_grid()：这就是所谓的预处理函数。 顾名思义，它设置存储三角测量的对象。 在后面的示例中，它还可以处理边界条件、几何形状等。
    \item setup\_system()：这是设置解决问题所需的所有其他数据结构的函数。 特别是，它将初始化 DoFHandler 对象并正确调整与线性代数有关的各种对象的大小。 该函数通常与上面的预处理函数分开，因为在时间相关的程序中，每当网格自适应细化时，至少每隔几个时间步调用它（我们将在步骤 6 中看到如何执行此操作）。 另一方面，在上面的预处理函数中设置网格本身仅在程序开始时完成一次，因此被分成自己的函数。
    \item assemble\_system()：这是计算矩阵和右侧内容的地方，如上面的介绍中详细讨论的。 由于对这个线性系统进行处理在概念上与计算其条目非常不同，因此我们将其与以下函数分开。
    \item solve()：这是我们计算线性系统 $AU=F$ 的解 $U$ 的函数。 在当前的程序中，这是一项简单的任务，因为矩阵非常简单，但只要问题不再那么微不足道，它将成为程序大小的重要组成部分（例如，参见步骤 20、步骤 22） ，或者在您对库有了更多了解后执行步骤 31）。
    \item output\_results()：最后，当您计算出解决方案时，您可能想用它做一些事情。 例如，您可能想要以可视化的格式输出它，或者您可能想要计算您感兴趣的数量：例如，热交换器中的热通量、机翼的空气摩擦系数、最大桥梁载荷，或者 只是一点的数值解的值。 因此，该函数是对解决方案进行后处理的地方。
\end{itemize}

所有这些都由单个公共函数（构造函数除外）（即 run() 函数）结合在一起。 它是从创建该类型的对象的位置调用的函数，并且是按正确顺序调用所有其他函数的函数。 将此操作封装到 run() 函数中，而不是从 main() 调用所有其他函数，可确保您可以更改此类中关注点分离的实现方式。 例如，如果其中一个函数变得太大，您可以将其分成两个，并且您唯一需要担心因此而改变的地方是在同一个类中，而不是其他任何地方。

如上所述，您将在以下许多教程程序中再次看到这种一般结构 - 有时函数名称的拼写有所不同，但本质上是按照功能分离的顺序。

\textbf{关于类型的注释}\\

deal.II 通过命名空间类型中的别名定义了许多整型。 （在上一句中，“积分”一词用作对应于名词“整数”的形容词。它不应与表示曲线或曲面下的面积或体积的名词“积分”混淆。 形容词“积分”在 C++ 世界中广泛用于“积分类型”、“积分常数”等上下文中。）特别是，在此程序中，您将在几个地方看到 types::global\_dof\_index： 整数类型 用于表示自由度的全局索引，即在三角剖分之上定义的 DoFHandler 对象内特定自由度的索引（与三角剖分顶部定义的特定自由度索引相对） 特定的细胞）。 对于当前的程序（以及几乎所有的教程程序），您将在全球范围内拥有几千到几百万个未知数（并且，对于 Q1 元素，您将在 2d 中的每个单元上本地有 4 个，在 3d 中的每个单元上有 8 个） ）。 因此，允许为全局 DoF 索引存储足够大的数字的数据类型是 unsigned int，因为它允许存储 0 到略大于 40 亿之间的数字（在大多数系统上，整数为 32 位）。 事实上，这就是 types::global\_dof\_index 的本质。

那么，为什么不立即使用 unsigned int 呢？ deal.II 在 7.3 版本之前一直这样做。 然而，deal.II 支持非常大的计算（通过步骤 40 中讨论的框架），当分布在几千个处理器上时，可能有超过 40 亿个未知数。 因此，在某些情况下，unsigned int 不够大，我们需要 64 位无符号整型。 为了实现这一点，我们引入了 types::global\_dof\_index ，默认情况下它被定义为简单的 unsigned int ，而如果需要，可以通过在配置期间传递特定标志将其定义为 unsigned long long int （请参阅自述文件）。

这涵盖了技术方面。 但还有一个文档目的：在库和基于其构建的代码中的任何地方，如果您看到使用数据类型 types::global\_dof\_index 的地方，您立即知道所引用的数量实际上是 全局自由度指数。 如果我们只使用 unsigned int （也可能是本地索引、边界指示符、材质 ID 等），则不会有明显的含义。 立即知道变量引用的内容也有助于避免错误：很明显，如果您看到类型为 types::global\_dof\_index 的对象被分配给类型为 types::subdomain\_id 的变量，即使它们都由 无符号整数，因此编译器不会抱怨。

从更实际的角度来说，这种类型的存在意味着在组装过程中，我们创建一个 $4\times4$ 矩阵（在 2d 中，使用 $Q_1$ 元素），其中包含我们当前所在单元的贡献，然后我们需要添加 该矩阵的元素到全局（系统）矩阵的适当元素。 为此，我们需要获取当前单元局部自由度的全局索引，为此我们将始终使用以下代码片段：
\begin{lstlisting}
cell->get_dof_indices (local_dof_indices);
\end{lstlisting}
其中 local\_dof\_indices 声明为
\begin{lstlisting}
std::vector<types::global_dof_index> local_dof_indices (fe.n_dofs_per_cell());
\end{lstlisting}
这个变量的名称可能有点用词不当——它代表“当前单元上局部定义的自由度的全局索引”——但保存此信息的变量在整个库中普遍以这种方式命名。

\section{step-3在计算的问题}

Step-3主要解决的是一个二维的Poisson方程。Poisson方程是一种常见的椭圆型偏微分方程，可以描述很多物理问题，例如电势、热传导等。

具体来说，示例程序通过有限元方法求解以下Poisson方程：

$\nabla^2u=f$

其中，$u$是未知函数（解），$f$是已知的右端项。

step-3的主要计算过程包括以下几个步骤：
\begin{enumerate}
    \item 创建网格：生成一个简单的二维正方形网格用于离散化问题域。
    \item 定义有限元空间：根据网格创建有限元空间对象，用于离散化解的空间。
    \item 定义边界条件：设置边界值条件，描述问题的边界约束。
    \item 构建线性方程系统：根据有限元离散化方法，构建线性方程的矩阵和右端向量。
    \item 求解线性方程：使用预处理共轭梯度（PCG）等求解方法，求解线性方程。
    \item 输出结果：将计算得到的解输出，以及其他感兴趣的信息。
\end{enumerate}

通过这些计算步骤，演示了如何使用DEAL.II库求解Poisson方程的基本过程，可以作为学习和理解有限元方法及其在实际问题求解中的应用的参考。

step-3中的具体参数如下：
\begin{enumerate}
    \item 网格参数：
    \begin{itemize}
        \item domain\_size: 设置问题域的大小，示例程序中使用的是2D正方形域，大小为$1\times1$。
        \item number\_of\_subdivisions: 设置网格的分割数，示例程序中使用的是$8\times8$。
    \end{itemize}
    \item 有限元空间参数：
    \begin{itemize}
        \item degree: 设置有限元的度数，示例程序中使用的是一阶连续有限元，即degree=1。
    \end{itemize}
    \item 边界条件参数：
    \begin{itemize}
        \item boundary\_value: 定义了边界上的数值，示例程序中使用的是DirichletBoundaryValues类，该类在特定的点上定义了边界值。
    \end{itemize}
    \item 线性方程系统参数：
    \begin{itemize}
        \item quadrature\_degree: 设置数值积分的度数，示例程序中使用的是2（二次精度）。
        \item system\_matrix\_sparsity\_pattern: 设置线性方程系统的稀疏模式，用于构建稀疏矩阵。
        \item solver\_control: 设置求解器的控制策略，例如迭代次数、收敛容限等。
    \end{itemize}
    \item 输出结果参数：
    \begin{itemize}
        \item output\_filename: 设置输出文件的名称，示例程序中使用的是solution.gpl。
        \item data\_out\_flags: 设置输出结果的内容，示例程序中输出解和梯度。
    \end{itemize}
\end{enumerate}

这些参数可以根据实际需求进行调整。例如，可以根据问题的复杂性和精度要求增加网格分割数和有限元的度数，以获得更准确的解。在实际应用中，还可以根据需要修改边界条件、添加其他方程项等。

\begin{figure}[h]
    \centering
    \includegraphics{step-3.pdf}
\end{figure}

\end{document}